
global data_in=20
global data_out=1180  
global bin=10
global binname=$bin 
global zoom_in_c300=80        
global zoom_out_c300=470
global range_c300   "incbin${binname}>=float(${zoom_in_c300}) & incbin${binname}<=float(${zoom_out_c300})"  

   
clear mata
matrix table=J(2,8,.)
global rr=1

use  "${density_data_dir}\density_${firm_name}_${name}_${year_name}_bin${binname}_in${data_in}_out${data_out}.dta", clear

global notch_s= ${cutoff}/${bin}  
global low_s=${low}/${bin}            //normalized; not include cutoff 		
global high_s=${high}/${bin}

global low_name    =-${low}
local  incbin_name ="incbin${binname}"
local  freq_name   ="numbin${binname}"   //original data

gen round50=0
replace round50=1  if mod(incbin${binname},50)==0

gen Z100=0
replace Z100=1     if mod(incbin${binname},100)==0   
	
	
		if "${firm}"=="foreign" {
	/*note: rounding at USD is included.*/	
		gen D10=0
		replace D10=1 if incbin${binname}==80 
		replace D10=1 if incbin${binname}==160 | incbin${binname}==170
		replace D10=1 if incbin${binname}==240 | incbin${binname}==250
		replace D10=1 if incbin${binname}==330 
		replace D10=1 if incbin${binname}==410 |  incbin${binname}==420
		
		gen D5=0
		replace D5=1 if incbin${binname}==120 
		replace D5=1 if incbin${binname}==280 | incbin${binname}==290
		replace D5=1 if incbin${binname}==370 
		replace D5=1 if incbin${binname}==450 |  incbin${binname}==460
		
     	gen D20=0
		replace D20=1 if incbin${binname}==160 | incbin${binname}==170
		replace D20=1 if incbin${binname}==330 
	
		gen D25=0
		replace D25=1 if incbin${binname}==200 | incbin${binname}==210
		replace D25=1 if incbin${binname}==410 | incbin${binname}==420	

		gen D50=0
		replace D50=1 if  incbin${binname}==410 |  incbin${binname}==420
		
		gen D_round=0
		replace D_round=1 if incbin${binname}==170
		replace D_round=1 if incbin${binname}==250
		replace D_round=1 if incbin${binname}==420
		replace D_round=1 if incbin${binname}==460
					
		gen D_c5=0
		replace D_c5=1 if incbin${binname}==370

		}		
		
		
			
keep if $range_c300


/* MATA ESTIMATION */		
mata:mata clear
mata:incbinname=st_local("incbin_name")   
mata:freqname=st_local("freq_name")    
mata:notch_s=strtoreal(st_global("notch_s"))  
mata:low=${low_s}
mata:high=${high_s}
mata:degree=${degree}
mata:nboot=${nboot}
mata:sf_s=${bin}
set seed 1234

* CALLING MATA PROGRAM TO ESTIMATE THE COUNTERFACTUAL, b, a*, z_u AND THE ELASTICITIES *
if "${density}" =="DP_q" & "${firm}"=="foreign"	{ 	
mata: notch_fie_ct_roundq(incbinname,freqname,notch_s,low,high,degree,nboot,sf_s)
}


/* POPULATING THE TABLE */
matrix table[1,1]=.
matrix table[1,2]=${cutoff}
matrix table[1,3]=${bin}
matrix table[1,4]=${degree}
matrix table[1,5]=${low}   //before normalize
matrix table[1,6]=${high}
matrix table[1,7]=${b}
matrix table[2,7]=${b_se}
matrix table[1,8]=${m}
matrix table[2,8]=${m_se}
	
noi di "b=" ${b}
noi di "b_se=" ${b_se}
noi di "m=" ${m}
noi di "m_se=" ${m_se}


/* GETTING THE BUNCHING PLOT */
svmat  y_actual_noroundsround 
svmat  y_counter
svmat  y_counter_noroundsround  

ren y_actual_noroundsround   numbin${binname}_noround          
ren y_counter_noroundsround  numbin${binname}_counter_noround  
ren y_counter                numbin${binname}_counter                 

save "${density_data_dir}\density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta", replace 
     
svmat table
keep table* 	
rename table1  Region
rename table2  cutoff
rename table3  bin
rename table4  degree 
rename table5  low
rename table6  high
rename table7  b
rename table8  m
	
tostring Region,replace
replace Region = " " if Region== "." 
replace Region="${name}"  in 1  
save "${density_data_dir}/Table1_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${bin}_in${zoom_in_c300}_out${zoom_out_c300}.dta",replace 
		

		
/*draw figs*/
global rr=1	
global b=strofreal(b[${rr}],"%9.2f")	
global b_se=strofreal(b[${rr}+1],"%9.2f")
global m=strofreal(m[${rr}],"%9.2f")	
global m_se=strofreal(m[${rr}+1],"%9.2f")
	
global graph_low=${low}+${cutoff}
global graph_high=${high}+${cutoff}


use "${density_data_dir}\density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta", clear



//figure 1A: exclusion region only
#delimit;
twoway 
(connected numbin${binname}                    incbin${binname}, sort clcolor(green) mcolor(green) msymbol(o) msize(large))
(line      numbin${binname}_counter          incbin${binname}, sort clcolor(maroon)  mcolor(maroon) msymbol(o) lwidth(medthick) lpattern(dash))
if incbin${binname}>=float(${graph_low}) & incbin${binname}<=float(${graph_high})
,
legend(label(1 "Observed")   label(2 "Counterfactual")  size(small))
xline(${cutoff}, lcolor(black) lwidth(med)) 
ytitle(Density, size(small)   margin(medium))  
xtitle(Registered Capital(unit:10,000 RMB), size(medsmall))xtitle(, alignment(top)) 						
xlabel(${graph_low} (${binname})  ${graph_high}, labsize(small))
yscale(r(0))	ylabel(#5) 
graphregion(fcolor(white)  ifcolor(white) color(white) icolor(white))
title("Density distribution of registered capital for", justification(center) color(black) size(medsmall)) 
subtitle( "${firm_name} firms in ${location} ${year_note}"  " ${cutoff_note}",  justification(center) color(black) size(med))
note("Note: bin=${bin}, degree=${degree}, b=${b}(${b_se}).",  size(small))
;
graph export 
"${density_graph_dir}/fitfig_density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_${low}_${high}_bin${binname}_${graph_low}_${graph_high}.pdf", replace;
#delimit cr










